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Abstract 

Particle-Mesh (PM) codes are still very useful tools for testing predictions of cosmological models in 
i cases when extra high resolution is not very important. We release for public use a cosmological PM N- 

0^ I body code. The code is very fast and simple. We provide a complete package of routines needed to set 

initial conditions, to run the code, and to analyze the results. The package allows you to simulate models 
O I with numerous combinations of parameters: open/flat /closed background, with or without the cosmological 

^ i constant, different values of the Hubble constant, with or without hot neutrinos, tilted or non-tiltcd initial 

t-H I spectra, different amount of baryons. 

■ Routines are included to measure the power spectrum and the density distribution function in your 

simulations, and a bound-density-maxima code for halo finding. We also provide results of test runs. A 
simulation with 256^ mesh and 128"^ particles can be done in a couple of days on a typical workstation (about 
70Mb of memory is needed) . To run simulations with 800'^ mesh and 256'^ particles one needs a computer 
with 1Gb memory and 1Gb disk space. The code has been successfully tested on an HP workstation and on 
a Sun workstation running Solaris, but we expect it should work on other systems. 
PsJ . The package can be downloaded from http:/ /astro. nmsu.edu/~aklypin/PM/PMcode.tar.gz A PostScript 

' version of this manual can be obtained from http:/ /astro. nmsu.edu/~aklypin/PM/PMcode.ps.gz 

. We provide this tool as a service to the astronomical community, but we cannot guarantee results. 

^ ■ 1 Introduction 

Q , There are many different numerical techniques to follow the evolution of a system of many particles. For 

^ ' earlier reviews see Hockney & Eastwood (1981) and Sellwood (1987). The most frequently used methods 

, for cosmological applications fall in three classes: Particle Mesh (PM) codes, Particle-Particle/Particle-Mesh 

' (P'^M) codes, and TREE codes. All methods have their advantages and disadvantages. 

PM codes use a mesh for the density and potential. As a result, their resolution is limited by the size of 
the mesh. The largest simulations done by the author have been on a 800'^ mesh with 3 x 256^ — 1.5 x 10^ 
particles. The SP2 parallel supercomputer at Cornell has been used to run simulations with a 1600^ — 
' 4.096 X 10^ mesh (Gross, 1997). There are two advantages of the method: i) it is fast (the smallest number 

of operations per particle per time step of all the other methods), ii) it typically uses very large number 
of particles. The latter can be crucial for some applications. There are a few variants of PM codes. A 
"plain-vanilla" PM was described by Hockney & Eastwood (1981), and this includes a Cloud-In-Cell density 
assignment and a 7-point discrete analog of the laplacian operator. Higher order approximations improve 
the accuracy on large distances, but degrade the resolution (e.g. Gelb (1992)). In an effort to reduce the 
order of approximation and to increase the resolution, Melott (1986) introduced the staggered mesh. It gives 
a better resolution on cell-size distances, but particles get self-forces (an isolated particle experiences a force 
from itself), which might be not a welcome feature. 

P'^M codes are described in detail in Hockney & Eastwood (1981) and in Efstathiou et al.(1985). They 
have two parts: a PM part, which takes care of large-scale forces, and a PP part, which adds a small-scale 
particle-particle contribution. The simulations usually have 64^-100^ particles. Because of strong clustering 
at late stages of evolution, the PP part becomes prohibitively expensive once large objects start to form 
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in large numbers. Significant speed is achieved in a modified version of tlie code which introduces subgrids 
(next levels of PM) in areas with high density (AP^M code of Couchman (1991)). With modification this 
code runs as fast as TREE code even for heavily clustered configurations (Couchman 1991). 

TREE codes are the most flexible code in the sense of the choice of boundary conditions (Appcl 1985, 
Barnes & Hut 1986, Hernquist 1987). They are also more expensive than PM: it takes 10-50 times more 
operations. Bouchet & Hernquist (1986) and Hernquist, Bouchet & Suto (1991) have extended the code for 
the periodical boundary conditions, which arc important for simulating large-scale fluctuations. 

Multigrid methods were introduced long ago, but only recently have they started to show a potential 
to produce real results (Anninos, Norman & Clarke 1994, Suisalu & Saar 1995, Kravtsov et al.l997). At 
present the most advanced and fastest multigrid code has been developed by Kravtsov et al.(1997). 

2 Equations and dimensionless variables 

The equations of motion of particles in expanding coordinates, which are used in our PM code, were presented 
by Kates et al.(1991). Different numerical effects (including resolution) were discussed in Klypin et al.(1996). 
We use comoving coordinates ofr particles x = x(t), which are related to proper coordinates by r = a{t)x, 
where a{t) = (1 + z)~^ is the expansion parameter. Instead of using peculiar velocity Vpec = ax we write 
the equations for particle momenta p: 

p = a^x, Vpec = p/a (1) 

This choice of "velocity" simplifies the equations of motion by removing a few terms with a/ a. It is also 
convenient to change the time variable from time t to the expansion parameter a. The equations governing 
the motion of particles are: 

dp \7(j) dx p 
da d ' da aa^ 

VV = 47rGf2^(t)aVcr(t)^ = 47rC?OoPcr,o-, 5 = ^Mll^, (3) 

a pb 



(2) 



d^/a = HofJi}o + ^^curv.OCl + flAfid^, ^0 + ^curv.O + ^A,0 = 1, (4) 

where Oo = ^m{z = 0), f^curv.o, and Qa.o sltc the density of the matter, effective density of the curvature, 
and the cosmological constant in units of the critical density aX z = 0. The curvature contribution is positive 
for negative curvature. 

Dimensionless variables (shown with tildas below) are defined by introducing the length of a cell of the 
grid xq and by measuring the time in units of I/Hq: 

x = a;ox, t = i/Ho, (5) 
Vpec = ixoHo)p/a, 4> = 4>{xoHQf (6) 



a3 87rG 



^0 (7) 
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Equations (2-5) can be rewritten in terms of dimensionless variables: 



F(„)V0, -=Fia)-^ (8) 



I^(P-I), (9) 
2 a 



X^2 

J7o + ilcurv,OCt + ^A.QO- 



where 

F{a) = Ho/a = (^ """ ' — ' j . (iQ) 

Equations (|| - ^) are solved numerically by the PM code. 

If L is the length of the computational box at z = 0, iVgrid is the number of grid cells in one direction, 
and A^row is the number of particles in one direction, which contribute a fraction Qq of the critical density, 
then the transformations from dimensionless variables given by the code to dimensional variables are given 

by 



v,..^i.oHof- =0.781^^ -l-J^^. (12) 
a s a A'gi.id/128 

Mass = 7Vparticie. -mi, mi = nopcr.o ( ^) = 1.32 ■ 105(r>o/i') f ) (13) 

3 Scheme of integration 

Equations (|| - ^ are solved using finite differences with a constant step in space Aa; — Ay — Az — 1 and 
a constant step in the expansion parameter Ao. We use the "leap-frog" scheme to advance coordinates and 
velocities from one moment to another. (In the following we drop tildas for all dimensionless variables.) At 
any moment a„ — ao + nAa, we have the coordinates x„ and the potential (/)„. Velocities Pn-1/2 a-re defined 
at a„_i/2 — an — Aa/2. The coordinates and the velocities for the next moment are found using: 



Pn+1/2 = Pn-1/2 - -P'(a«)V(/)„Aa, 
, F{a„+i/2) 

x„+i = x„ H 2 Pn+i/2Aa (14) 

^n+1/2 

In order to solve eq.(^) we approximate the Laplacian operator using the 7-point "crest" template: 

V^^ « 0i±lj\fc + 4>i,j±lM + 0jj,fc±l - Q4>i,j,k, (15) 

where (z, j, k) — 1, iVgiid. This leads to a large system of linear equations relating unknown variables (pij^k 
with known right-hand side of the discrete form of the Poisson equation 3flo{pi,j,k ~ l)/2a. The system of 
equations is solved exactly by the FFT technique. 
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The density on the mesh Pij^k is obtained from particle positions using the Cloud-In-CcU method. In 
order to find the "acceleration" g = — V(/> for each particle, the gravitational potential is differentiated on 
the mesh: 

Qx = -(0i+i,j,fe - 0i-i,j,fe)/2, Qy = Qz = ••• (16) 

Then the acceleration is interpolated to the position of the particle using a three-linear interpolation. This 
scheme for the force interpolation (the same interpolation as in the density assignment) is very important 
because it does not produce a force acting on the particle itself. (Thus, an isolated point does not produce 
a force at the position of the particle). While this might sound like a natural condition for any realistic 
method, only two methods - PM and TREE - do not have this self-force. In P"^M the effect is minimized. In 
the case of multigrid methods the self-force cannot be excluded - only minimized. Typically this is achieved 
by placing extended buffer zones around regions with high resolution (e.g. Kravtsov 1997). No precautions 
were made in AP'^M method, which might result in spurious effects in regions were multi-level grids are 
introduced. 

Thus, the main scheme of the PM method consists of the following four blocks repeated every time step: 

• Find density on the mesh using the Cloud-In-Cell technique. 

• Solve the Poisson equation using two three-dimensional FFTs. 

• Advance velocities and coordinates of the particles. 

• Advance time and print results. 

4 Format of data 

In order to have the best possible resolution, most of the available computer memory is allocated to the 
density/potential grid. The Poisson solver is organized in such a way that only one large mesh is needed. 
Particle coordinates and velocities are kept on disk and are read into memory in large portions when necessary. 
This reading/writing of particles results in a small overhead - typically 5-10% of the total cpu time. If this 
overhead is an issue, the code can be easily adjusted to keep all particles in memory. This is always the 
case for a parallel version of the code. Particles are divided into "species" with constant mass of a particle 
for each species. Each species is kept in a separate file. Information which describes the run (such as the 
number of particles, omegas, and current time) is written in a separate header file. 

Each file with particle data is a FORTRAN direct-access file with the number of records equal to 
the number of particles in one direction A^row Each record has coordinates and velocities for a "page" 
of particles A^pago = The "page" of particles is read into a common block, which has the struc- 

ture: A:(A^page),l'(A^page),-Z'(Afpage),14(A^page), Vy(Arpage),14(A'page) The particle files and the header file 
are needed for continuation of the run or for the data analysis. The following diagram shows the structure 
and names of the data files: 
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C3CRD.DAT Header Text-of-Header, a, ainit, Aa, Step,. . . 



C3crsO.DAT 


Set 


Page 1 


Xi, 


,X2, ■ 


• ' ^Npage ) 2/1 ?• • 


■ Zi,.. 




• ^Npage 






Page 2 


Xi, 


:X2, ■ 


• • ^Npage ) 2/1 ?• • 


■ Zi,.. 




• ^Npage 






Page N^ow 


Xi, 


:X2, ■ 


• • ^Npage ) 2/1 ?• • 


■ Zi,.. 




• ^Npage 


C3crsl.DAT 


Set 1 


Page 1 


Xi, 


.X2, ■ 


• • ^Npage ) 2/1 ?• • 


■ Zi,.. 


.. Kl,.. 


• ^Npage 






Page 2 


Xi, 


.X2, ■ 


• • -^Npagc 7 Ul 1 • • 


■ Zi... 


,. I4i,.. 


• ^Npagc 






Page iVrow 


Xi, 


.X2, ■ 


• • •^Npage ) Uli • • 


■ Zi,.. 




■ ^Npage 



File Name Description Content of the file 

Not a part of the file 

Thus, the memory required to run the code is about N^^^^ + QN^ow memory words or 64Mb(7Vgrid/256)^ 
+0.375Mb(-/Vrow/128)'^ if single precision arithmetic is used. The total amount of disk space is 
48Mb(7Vrow/128)^ per each set of "species". 

5 Initial conditions: CDM and CHDM models 

We use the Zeldovich approximation to set initial conditions. The approximation is valid in mildly nonlinear 
regime and is much superior to the linear approximation. Wc slightly rewrite the original version of the 
approximation to incorporate cases (like CHDM) when the growth rates b{t) depend on the wavelength 
of the perturbation |fc|. In the Zeldovich approximation the comoving and the lagrangian coordinates are 
related in the following way: 



x = q-a^6|fe|(f)Sjfcj(q), p = -aa^ ^ 6|ft|(t) | ^ j Sj^j (q), (17) 



where the displacement vector S is related to the velocity potential $ and the power spectrum of fluctuations 

P{\k\): 

S|fe| (q) = Vg$|fc| (q), = X! cos(kq) + 6k sin(kq), (18) 

k 

where a and b are gaussian random numbers with the mean zero and dispersion cr^ = P{k)/k^: 

«.= V«-^=f^, ^. = VpM-^=P^. (19) 

The parameter a, together with the power spectrum P{k), define the normalization of the fluctuations. 

We estimate the power spectrum P{k) for a wide range of cosmological models using a Boltzman code 
(Holtzman 1989). As compared with the original version of the code, the current version allows for more 
accurate estimates at high wavenumbers. For each cosmological model the numerical data points were fitted 
using the following fitting formula: 
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P(k) = k^eMPi) .201 

^ ' (l+P2fcl/2+P3fc + P4fc3/2^_Pgfc2)2Pe- ^ ' 

The coefficients Pi are presented in the file cdm.fit for a variety of models. The errors of the fits are smaller 
than 5% in the power spectrum. The top panel in Figure 1 shows the errors of the fits for CDM models 
{Qq — 1) with a Hubble constant H = 50km/s/Mpc. Errors at a level of ^ 2% level at fc ~ 3/i Mpc~^ 
and at fc ~ 30/i Mpc""'^ are due to small mismatch in approximations used at high wavenumbers. The fits 
smooth out the jumps and, thus, provide better approximations to the real power spectra at those large 
wavenumbers. The waves around k ~ O.lh Mpc~^ are due to acoustic oscillations in baryons. They are 
larger for larger flb/fla ratios. For very small flb/flo the errors introduced by using the fits are extremely 
small. Thus, if one can neglect (or smooth out) the acoustic oscillations, the maximum errors of our fits 
are expected to be smaller than 1-2% in the power. The comparison of some of our power spectra with 
the results from COSMICS (Bertschinger 1996) support our conclusion. We recommend the use of the fits 
whenever it is possible. 

The power spectrum of cosmological models is often approximated using a fitting formula given by 
Bardeen et al.(1986, BBKS): 

P{k) - fc"r2(fc), T(fc) = Mi±^:^[i + 3.89g + {W.lq? + (5.4q)3 + (6.71(7)4]-i/4^ (21) 

2.34g 

where q — k/{^Qh? Mpc~^). Unfortunately, the accuracy of this approximation is not great. Peacock & 
Dodds (1994) modified the fit using another relation between q and k: 

q = k/{noh^ exp(-217b) Mpc"^). (22) 
This approximation was criticized by Sugiyama (1995), who introduced a better scaling for low-f2o cases: 

kiTcMB/2.7Kf 



Qoh'^ cxp{-Qb - v/VO^^b/^o) Mpc~ 



(23) 



These approximations have been frequently used in a number of publications (e.g. Liddle et al.l996). The 
bottom panel of Figure 1 shows the ratio of the power spectrum given by this approximation to the power 
spectrum obtained from our fits for several choices of baryon fraction. For comparison, we also present the 
error of the eqs.(|T[^) relative to the power spectrum obtained by COSMICS for fif, ~ 0.05 (triangles), 
showing the good agreement of our results with those of COSMICS. In all cases, there is a large decline 
(around 20% in power) between a peak at fc = 0.2h Mpc~^ and small scales k ~ (10 — 30)h Mpc~^. This 
decline was noticed by Hu & Sugiyama (1996), who studied the small-scale perturbations. Note that if we 
take TcMB = 2.70K instead of 2.726K, than the peak of the error at A: = 0.2 increases up to 15%. The error 
in the power is rather small for small k < O.lh Mpc~^ and for a realistic amount of baryons fib ^ 0.07. One 
can easily miss it if instead of an error of the power spectrum, one plots the transfer function in a double 
logarithmic scale. But the error is very significant for galactic-scale events. It can result in serious errors in 
the epoch of galaxy formation or in the amount of gas in damped Ly-a clouds at high redshifts. 

Hu & Sugiyama (1996) recommend changing the last parameter in the BBKS fit from 6.71 to 6.07. 
We do not find that this correction gives an accurate fit to our spectrum. We find that the following 
approximation, which is a combination of a slightly modified BBKS fit and the Hu & Sugiyama (1996) 
scaling with the amount of baryons, provides errors in the power spectrum smaller than 5% for the range of 
wavenumbers k = (10^** - 40)/i Mpc~^ and for ilb/fio < 0.1: 
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Figure 1: (Top) Errors of the fits eq.(|20|) for the CDM models {Qq = 1) with a Hubble constant H = 
50km/s/Mpc. Errors at the ~ 2% level at fc ~ 3/i Mpc~^ and at fc ~ 30h Mpc~^ are due to a small mismatch 
in approximations used at high wavenumbers. The fits smooth out the jumps and, thus, provide better 
approximations to the real power spectra at these large wavenumbers. The waves around fc ~ O.l/i Mpc~^ 
are due to acoustic oscillations in baryons. (Bottom) The differences between the power spectrum given by 
the BBKS approximation and the power spectrum obtained from our fits. Triangles show results obtained 
using COSMICS for fib ^ 0.05 
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P{k) = fc"T2(fc), 

T{k) = ^I^il±^pl[l + 13q+{10Mf + {WAqf + {6.51q)T'/\ 

ai = (46.90o/i')°'™[l + (32.1f^o/i')-°-'^'], a2 = (12fio/i')°-^'^[l + (45fio/i')-°'''] (24) 
Figures 2 and 3 show errors of the fits for the CDM and for the ACDM models. 



6 Finding Halos with Bound Density Maxima code 

Finding halos in dense environments is a challenge. The most widely used halo-finding algorithms - the 

fricnds-of-fricnds (e.g., Efstathiou et al.l985) and the spherical overdensity algorithm (e.g., Laccy & Cole 
1994; Klypin 1996) - are not acceptable (Gelb & Bertschinger 1994, Summers et al.l995). The friends- 
of-friends (FOF) algorithm merges together apparently distinct halos if the linking radius is too large or 
misses some of the halos if the radius is small. Adaptive FOF (van Kampcn 1995) seems to work better. 
But we find that it is difficult to find an optimal scaling of the linking radius with the density. We have 
developed a related algorithm (Klypin, Gottlober, Kravtsov 1997), which we call "hierarchical friends-of- 
fricnds" . Because it uses all linking radii, it docs not have the problem that the adaptive FOF algorithm 
has. The algorithms, either adaptive or hierarchical, can not work because they pick up many fake halos in 
very dense environments. Klypin et al.(1997) supplement the hierarchical FOF algorithm with an algorithm 
which checks if halos existed at previous moments. The algorithm which finds halos as maxima of mass 
inside spheres of a given overdensity works better than plain FOF, but no fixed overdensity limit can find 
halos in both low and high density environment. The DENMAX algorithm (Gelb & Bertschinger 1994) 
and its offspring, SKID (Governato et al.l997), make significant progress - they remove unbound particles, 
which is important for halos in groups and clusters. Recently, Summers et al.(1995) tried to perfect the idea 
of Couchman & Carlberg (1992) to trace the history of halo merging and to use it for halo identification. 
Starting at an early epoch. Summers et al. identify halos using the FOF algorithm with a linking radius 
corresponding to the "virial overdensity" of 200 and then trace particles belonging to halos at later times. 
It appears that it is impossible to make a working algorithm because halos interact too violently. A large 
fraction of mass is tidally stripped from some halos and a large fraction of mass is accreted. Some of the 
problems that any halo finding algorithm faces are not numerical. They exist in the real Universe. We select 
a few typical difficult situations. 

1. A large galaxy with a small satellite. Examples: LMC and the Milky Way or the M51 system. 
Assuming that the satellite is bound, do we have to include the mass of the satellite in the mass of the large 
galaxy? If we do, then we count the mass of the satellite twice: once when we find the satellite and then 
when we find the large galaxy. This does not seem reasonable. If we do not include the satellite, then the 
mass of the large galaxy is underestimated. For example, the binding energy of a particle at the distance of 
the satellite will be wrong. The problem arises when we try to assign particles to different halos in an effort 
to find masses of halos. This is very difficult to do for particles moving between halos. Even if a particle 
at some moment has negative energy relative to one of the halos, it is not guaranteed that it belongs to the 
halo. The gravitational potential changes with time, and the particle may end up falling onto another halo. 
This is not just a precaution. This actually was found very often in real halos when we compared contents 
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Figure 2: Errors of the 
amount of baryons. 
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Figure 3: The same as Figure 2, but for the ACDM models. 
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of halos at different rcdshifts. Interacting halos exchange mass and lose mass. We try to avoid tlie situation: 
instead of assigning mass to halos, we find the maximum of the "rotational velocity" , sjGMjR, which is 
observationally a more meaningful quantity. 

2. A satellite of a large galaxy. The previous situation is now viewed from a different angle. How can we 
estimate the mass or the rotational velocity of the satellite? The formal virial radius of the satellite is large: 
the big galaxy is within the radius. The rotational velocity may rise all the way to the center of the large 
galaxy. In order to find the outer radius of the satellite, we analyze the density profile. At small distances 
from the center of the satellite the density steeply declines, but then it flattens out and may even increase. 
This means that we reached the outer border of the satellite. We use the radius at which the density starts 
to flatten out as the first approximation for the radius of the halo. This approximation can be improved by 
removing unbound particles and checking the steepness of the density profile in the outer part. 

3. Tidal stripping. Peripheral parts of galaxies, responsible for extended flat rotation curves outside of 
clusters, are very likely tidally stripped and lost when the galaxies fall into a cluster. The same happens 
with halos: a large fraction of halo mass may be lost due to stripping in dense cluster environments. Thus, 
if an algorithm finds that 90% of mass of a halo identified at early epoch is lost, it does not mean that the 
halo was destroyed. This is not a numerical effect and is not due to "lack of physics" . This is a normal 
situation. What is left of the halo, given that it still has a large enough mass and radius, is a "galaxy" . 

We have developed our halo-finding algorithm (Klypin et al.l997) having in mind all these problems. The 
boimd-density- maxima (BDM) algorithm stems from the DENMAX algorithm (Gelb & Bertschinger 1994). 
Just as in DENMAX, the BDM algorithm first finds positions of the density maxima on some scale and 
then removes unbound particles inside the halo radius. However, the algorithm finds maxima and removes 
unbound particles in a way different from DENMAX. The algorithm can work by itself or in conjunction with 
the hierarchical FOE. In the latter case, it takes positions of halos from the hierarchical friends-of-friends, 
and then removes unbound particles and finds parameters of halos. 

In order to find positions of halos we choose a radius r^p of a sphere for which we find maxima of mass. 
This defines the scale of objects we are looking for, but not exact radii or masses of halos. The radius of 
a halo can be either larger or smaller than r^p, but distances between halos cannot be smaller than r^p. 
We place a large number of spheres in the simulation box. The number of the spheres is typically an order 
of magnitude or more larger than the expected number of halos. For each sphere we find the mass inside 
the sphere and the center of the mass. The center of the sphere is displaced to the center of mass and the 
new mass and the center of mass is found. The process is iterated until convergence. Depending on specific 
parameters of the simulations, the number of iterations ranges from 10 to 100. This process finds local 
maxima of mass within sphere of radius rgp. 

The efficiency of finding local maxima of mass depends on how the spheres are chosen. In the present 
version of the code two algorithms were implemented. (1) A small fraction of all particles are chosen as 
centers of the spheres. The code asks you to enter Ngeed, the "Number of particles for initial seeds". Then 
it will select every Nparticies / N seed particle as a center of a sphere. (2) Additional spheres will be placed in 
regions with relatively low density. The whole simulation box is divided into a mesh of large cells. The size 
of the cells is defined by the variable "Cell" in PMhalos.f. The "Cell" is typically equal to one or two PM 
cells. If a cell has many particles in it ("neighbors"), than some of them (not more than 3) will be chosen 
as centers of spheres. The code asks you to enter the minimum number of neighbors: "Number of neighbors 
for a seed" . 

In some cases one would need to improve the location of the halo. An example is if one is looking for 
groups of "galaxies" but also would like to have the groups always centered on a galaxy. The search radius 
for the groups may be chosen to be, say 500 kpc. Additional iterations with a smaller radius of the sphere 
will find the galaxy-size halo closest to the center of mass of the group and place the center of the group at 
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the "galaxy" . In the BDM code this option is rcahzcd in the fohowing way. The code asks you to enter the 
"smaller radius for final hales". If this radius r small is not equal to the search radius r^p, the code will do 
additional iterations by gradually changing the search radius from r^p to r small- If f small = fsp, no additional 
iterations are made. 

Some of the density maxima will be found many times because in the process of maximizing the mass 
some of spheres converge on the same local maximum. Spheres which find the same maximum are called 
"duplicates" . We remove duplicates and keep only one halo for each maximum. Halos with too small 
number of particles (typically 5-10) and halos with too low central overdensity are removed from the final 
list. Parameters which control the removal are supplied by the user. 

Once centers of potential halos are found, we start the procedure of removing unbound particles and 
finding the structure of halos. We place concentric spherical shells around each center. For each shell we find 
the mass of the dark matter particles, the mean velocity, the velocity dispersion relative to the mean, and 
the maximum of the rotational velocity VJnax = ^/ GM{r) jr \rna:r- In order to determine whether a particle 
is bound or not, we estimate the escape velocity at the distance r of the particle from the halo center: 



KUeW « (2.15 X y_)2Mld_^!V!W^, (25) 

V' / ' max] 



where Vmax is the radius of the maximum of the rotational velocity. This expression for the escape velocity 
is valid for a halo with the Navarro-Frcnk- White density profile. If the velocity of a particle is larger than 
the escape velocity, it is assumed to be unbound. We estimate the maximum rotational velocity Vmax and 
radius of the maximum Tmax using the density profile for the halo. Because Vmax and Vmax must be found 
before the unbound particles arc removed and because the mean velocity is also found using all particles 
(bound and unbound), the whole procedure cannot be done in one step. We start by artificially increasing 
the value of the escape velocity by a factor of three. Only particles above the limit are removed. We find a 
new density profile, new mean velocities, and The escape velocity is again increased, 

but this time by a smaller factor. The procedure is repeated 6 times. The last iteration does not have any 
extra factors for the escape velocity: all unbound particles are removed. Examples of halos identified by the 
code are presented in Klypin ct al.(1997). 

Finding a halo radius is straightforward for isolated halos: increase the radius of sphere untill the over- 
density inside the sphere is equal to some limiting value provided by the user. For halos inside groups or 
clusters (halos inside halos) the algorithm is more complicated. It consists of three steps: (i) It starts with 
finding the radius of given overdensity limit (as for an isolated halo), (ii) Then, the algorithm checks how 
the mean overdensity inside given radius changes with the radius. It starts going from the very central 
shell outwards. If the mean overdensity stops declining, the algorithm assigns the radius to the halo radius, 
(iii) Now the algorithm goes from this radius inwards and checks the slope n of the overdensity profile: 
{M{r)/{M{r))) oc r~". If the slope is too shallow (mass increases too fast with the radius), the radius of 
the halo is decreased because most of the mass at this radius does not belong to the halo. The radius is 
decreased untill the slope is steep enough. The limit for the slope is equal to n = 1. This is a rather mild 
slope: for an isothermal sphere one expects n = 2; for the Navarro-Frenk- White profile n = 1.7 — 2.5. 

The BDM code will ask you to enter many parameters. A typical dialog may look as follows: 

Enter Min. Center Overdensity for Halos => 340. (1) 

Enter Overdensity Threshold for Halos => 340. (2) 

Enter Minimum halo mass in Msim/h 5.e-|-9 (4) 

Enter comoving search radius(Mpc/li) => 0.050 (5) 

Enter smaller radius(Mpc/h) of final halos=> 0.030 (6) 
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Enter min. radius for halos(Mpc/h) => 0.030 (7) 

Enter fraction of DM particles (1/4,1/2, !)=> 1 (8) 

Enter rejection velocity limit (V/Vescape)=> 1.0 (9) 

Distance to check for Velocity duplicates => 0.300 (10) 

Define duplicates if (vl-v2)/Vrms < 0.10 (11) 

Enter Comoving Box size(Mpc/h) => 20. (12) 



In lines (1-2), enter the minimum overdensity for halos. If both values are equal, the code will find halos 
with the overdensity above the limits you provided. You may choose only those halos which have higher 
central density, if you enter a larger number in the first line. Only halos with mass larger than the value 
entered in the third line and radius in the line 7 will be kept. For debugging of the code, you may read 
only 0.25, or 0.5 of all particles (line 8). If you would like to remove unbound particles as described above, 
enter 1 in line 9. In order to ignore the option (no removing of unbound particles) enter a number larger 
than 5. Lines 10 and 11 are used for two parameters needed for additional screening of duplicates. In case 
of very large halos, when the density profile in the center of the halo is rather flat, the iteration of spheres 
stops when the spheres arc still far one from the other. The spheres have found the same object: masses, 
radii, velocities of the " halos" are very close, but their positions are slightly different. Because the fake 
halos have very close velocities, they can be identified. The parameter in the line 10 defines the maximum 
distance (in units of h~^Mpc) within which the code will look for the duplicates. If the difference in velocities 
V {vxi — 'Vx'i)^ + (-J/) + (-z) is less than the parameter in line (11) as measured in units of maximum of the 
rms velocities for particles in the halos, only largest will be kept. You may ignore the option by entering two 
zeros. 

7 How to compile and run the code 

The code consists of several different FORTRAN programs: 

• PM_to_ASCII.f 

• PMhalos.f 

• PMmain.f 

• PMmodelCHDM.f 

• PMmodels.f 

• PMpower.f 

• PMselect.f 

• PMstartCDM.f 

• PMstartCHDM.f 

A Makefile is provided to allow easy compilation of all of the programs. You should edit the Makefile 

and put in your preferred compilation flags for: 

• Optimization. On both our Sun machines running Solaris, and our HP running HPUNIX, we use -03. 



13 



• The routine PMmodels.f should be compiled using double precision for all "real" variables and 
constants. On the Sun, this can be accomplished using the -R8 flag, on the HP with -dbl. 

The general scheme for running the code and analyzing results is as follows: 

1. Set initial conditions using PMmodels and PMstartCDM. 

2. Run the PM code using PMmain. 

3. Analyze the results using PMpower, and PMhalos. PM_to_ASCII will scale your results to "nor- 
mal" units. 

More details on each step are provided in subsequent sections. 
Some important parameters and variables which are used in the codes, for which you will either be prompted 
or may wish to change before compilation. Please note the last item which may require you to make 
a change in the code! 

• AEXPN = expansion parameter (= 1/(1 + z)). AEXPO is the expansion parameter at initial moment. 

• NROW = number of particles in one dimension (= 2"). The total number of particles is equal to 
NROW^. Particles are stored in direct-access files. Each record of the files contains NROW^ particles. 

• NGRID = size of the computational mesh (= 2™). The total number of cells is NGRID^. If you 
need to change either NROW or NGRID, the only place where you need to make changes is the file 
PMparameters.h. 

• ASTEP = step in the expansion parameter a. Time integration is done with a constant step in the 

expansion parameter: da =ASTEP. 

• ISTEP = current integration step 

• Nspecies = number of additional species of particles. This is curretly used for the CHDM models. For 
plain CDM, ODM, or ACDM models Nspecies = 0. 

• OmO,OmlO,Ocurv = densities of matter, cosmological constant, and curvature at ^ = 0. 

• bubble = the Hubble constant in units of lOOkm/s/Mpc 

• NACCES = length of a record of direct-access files with coordinates and velocities of particles (e.g., 
PMcrsO.DAT). For some computers the length of the record is counted in bytes, for some it is in 
machine words. If you get errors when running the code that indicate problems with accessing the 
data files, multiply NACCES by 4 in file PMauxiliary.f, routine RDTAPE, and in files PMstart...f. 
If the length of the data files is too long (it should be 24 x NROW^ bytes), divide NACCES by 4. 



8 How to set initial conditions 

Two programs need to be run in order to set initial conditions. If the model you are interested in do not have 
massive neutrinos, the programs are PMmodels and PMstartCDM. For models with massive neutrinos, 
the corresponding programs are PMmodelCHDM and PMstartCHDM. Additionally, if you wish to run 
the PM code for a model for which the initial power spectrum is not included in the package, it is possible 
for you (with a bit more work) to specify whatever power spectrum you want. 
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8.1 Models without massive neutrinos 



A range of initial spectra which have been computed by our Boltzmann code and fit with our fitting function 
(ea.^O|) are available. The fitting parameters are tabulated in the file cdm.fit. You will need to look in this 
file and determine the line number of the model which you are interested in. The first program (PMmodels) 
will ask to give you the following parameters: 

• Name of a file where it writes parameters of the model and normalized power spectra. 

• (Tg = rms density perturbation in a sphere of radius 8/i^^Mpc and the slope of the power spectrum at 
long 

• Line number in the file cdm.fit 

• Size of the simulation box, number of particles in ID, and the redshift at which you would like to start 
your code (or to get the power spectrum). 

• The program will ask you if you need to have a file with all input parameters needed to set initial 
conditions. If you need the file, answer "Yes", and it will create file InStart.dat, which you will give 
as input to the second code PMstartCDM. It will also ask you few questions needed to produce the 
file. If you answer "No" , the code finishes its work by creating the file with the name you gave it in the 
first line. The file has parameters of the model (all fi's, the Hubble constant, the age of the universe, 
growth rate of density 5, and d{\n 5) / d{ln a) at different redshifts. It also gives the bulk velocity of a 
sphere of radius 50h~^yipc and normalized power spectra of dark matter at redshift zero (a = 1) and 
at the redshift you provided. 

• If you answered "Yes" for the previous question (you needed to have an input file for PMstartCDM), 
you will be asked to provide the following information: 

— A string of up to 45 characters ("header"). The header is not used for simulations, but it is useful 
to label your simulation. You can provide any information you want. This header will stick to 
your run. All files with snapshots of your simulation and all files with analysis of your simulation 
will have the header. Experience of running many simulations shows that one never has enough 
information describing details of a simulation done some time ago to identify that simulation later. 
Use the header to identify your run. 

— Step in the expansion parameter da. This defines how many integration steps the code will do 
untill it runs to the redshift zero. If you would like to make N steps and you start at redshift 
z, da = (1 — The step da should be (significantly) smaller than the initial expansion 
parameter Oinit = -j^- 

— Seed for random numbers. Use any integer number in the range 1 — (2'^^ — 1) (2^^ « 2.1478 x 10^). 

Next you will run PMstartCDM which is the program that actually generates initial conditions for the 
PM code. It will generate two files, which PMmain reads and updates: PMcrd.DAT (information on the 
cosmological model and parameters of the run) and PMcrsO.DAT (coordinates and velocities of particles). 
If you run a CHDM simulation, you will have more data files with data for hot neutrinos. PMstartCDM 
will ask you to provide some parameters. If you created file InStart.dat using the program PMmodels, 
simply provide the file as input: PMstartCDM < InStart.dat. 
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8.2 Models with massive neutrinos 

In order to set initial conditions for a CHDM model, run PMmodelCHDM and then PMstartCHDM. 
The codes will ask you similar questions as for non-CHDM models. Because setting initial conditions for the 
CHDM model is more complicated, we provide only one particular CHDM variant: a model with two equal 
mass neutrinos with total contribution il^ = 0.20, h = 0.5. Initial conditions can be set only at redshift 
z = 30. 

8.3 Initial conditions for arbitrary initial power spectra 

There are two ways of building initial conditions for models which arc not provided with the package. 

1. You may add a line to cdm.fit with fitting parameters for your model and with parameters of ap- 
proximation for the power spectrum. This requires that your model is well fit by our fitting formula 
(eq. ^0|). The format of the cdm.fit file is described in the file. See also the routine TRUNF(k) for 
details of the approximation of the power spectrum of perturbations. 

2. You may change the functions TRUNF in PMstartCDM.f and Pk in PMmodels.f to return what- 
ever initial spectra you desire. 

9 How to run the PM code 

After you generate the initial conditions, you can start running the code. Check if two files PMcrd.DAT 
and PMcrsO.DAT were generated and are in your directory: PMmain will read the files. It will also 
overwrite the files when it finishes. So, if you need to have the original files, please make copies before you 
start running PMmain. 

10 How to get coordinates and velocities of particles 

To get the final coordinates and velocities, run PM_to_ASCII. This will convert the input files PM- 
crd.DAT and PMcrsO.DAT into readable output files. You will specify the output file name. 

11 How to get power spectrum and density distribution 

You can get the power spectrum and density distribution of your output by running PMpower. This will 
read the raw output files from PMmain, PMcrd.DAT and PMcrsO.DAT. The output will go into a file 
called Spectrum.DAT. 

12 How to find halos using the Bound-Density-Maxima code 

Run PMhalos. The code will ask you many questions. It will produce two files with the final results. 
"Catalog.DAT" contains detailed information about all halos found by the code. After a rather long preamble, 
data on the halos follow. Each halo has a "header", which gives global parameters: coordinates, velocities, 
mass, and so on for the halo (format is given in the preamble). "Catshort.DAT" has a shorter preamble and 
has only a list of the halo headers. 



16 



Tips: The code may miss some halos if one chooses wrong parameters! 

1. The number of spheres (seeds) should be very large: 100,000 - 150,000. If it is too small, the code 
misses some of small halos (but not the big ones). 

2. The mimbcr of particles in each shell for the halo profile should be large: 5-6 per shell. The code needs 
the density profile of a halo to find radius, escape velocity and so on. It may get confused if the profile 
is too noisy. 

3. Radius of the first bin for the profile should not be too small: not less than 1/2 of your force resolution. 
It should also contain few particles (> 2). 

4. Outer radius should be large. If density docs not decline enough (flat profile), the code thinks that 
this halo is a fluke. If your radius is too small such that only the central part of a potential halo fits 
in the radius, and the density gradient is not steep, the code will kill your halo. 

13 Examples 

We make available two examples of runs of the PM code with data analysis which you can use for comparison 
with your own results if that is desired. The first example has files for the test with 32*^ particles and 128^ 
mesh. The second example is for 128^ particles and 256^ mesh. Both tests were done for a ACDM model 
with h = 0.7. 

These are available as two separate tar files with all results or you can also grab individual result files. 
Since the output data files for the second example are quite large, the tar file for this case does not include 
these data files. The full set of files from the 32'' example can be downloaded from http:/ /astro. nmsu.edu/ 
~aklypin/PAI/TEST32xl28.tar.gz. The files from the 128^ example (without data files) can be downloaded 
from: http://astro.nmsu.edu/-aklypin/PM/TEST128x256.tar.gz. 

Individual output files from the two cases can be downloaded from: http://astro.nmsu.edu/~aklypin/ 
PM/TEST32xl28 and http://astro.nmsu.edu/~aklypin/PM/TEST128x256. We recommend looking at the 
output plots (*.ps.gz). 
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